Non-perturbative theoretical description of two atoms in an optical lattice with 

time-dependent perturbations 
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A theoretical approach for a non-perturbative dynamical description of two interacting atoms 
in an optical lattice potential is introduced. The approach builds upon the stationary eigenstates 
found by a procedure described in Grishkevich et al. [Phys. Rev. A 84, 062710 (2011)]. It allows 
presently to treat any time-dependent external perturbation of the lattice potential up to quadratic 
order. Example calculations of the experimentally relevant cases of an acceleration of the lattice 
and the turning-on of an additional harmonic confinement are presented. 



I. INTRODUCTION 

Triggered by the creation of the first Bose-Einstein con- 
densates [H, [1| , the field of ultracold atoms has experi- 
enced many major advancements. Nowadays it is not 
only possible to steer and observe many-body effects like 
the Mott-insulator superfluid phase transition but 
also to manipulate single atoms in an optical lattice (OL) 
or a dipole trap d, 0] • 

A key technology is the dynamical variation of the 
trapping potential that allows, e.g., for a cooling of the 
system by transferring hot atoms to non-trapped contin- 
uum states [|[ . In a recent work by some of us it has been 
proposed how to perform quantum computations in an 
OL just by manipulating the depth of single lattice sites 
and by shaking the optical lattice to drive transitions 
between different Bloch bands @ . For a full understand- 
ing of the unde rlying dynamical processes of any multi- 
band system p| [lol - [l3 | the application of the usually em- 
ployed single-band Hubbard model is insufficient. Here, 
a numerical approach is presented, that solves the full 
time-dependent Schrodinger equation of two interacting 
atoms in a single-well or multiple-well lattice, which can 
be perturbed by any additional time-dependent potential 
up to quadratic order. While the types of perturbations 
can be easily extended, the currently implemented types 
already allow for studying many experimentally relevant 
situations. For example, an acceleration of an OL or a 
periodic driving as realized in [ill [l4| results in a lin- 
ear perturbation of the lattice. The m anip ulation of the 
barrier hight between two lattice sites [10[ or a variation 
of the global confinement, e.g. by a MOT [l5j . can be 
simulated by adding a harmonic perturbation. 

The general problem of a precise description of inter- 
acting atoms in trapping potentials is the existence of two 
very distinct length scales: that of the short-range Born- 
Oppenhcimer interaction (some 100 a.u.) and that of the 
trapping potential (some 10 000 a.u.). The employed ba- 
sis functions have to cover the highly oscillating behavior 
in the interaction range and the slow variation due to the 
trap. The use of an uncorrelated basis such as a regular 
grid or products of single-particle solutions is therefore 
impractical. A method to avoid the length scale prob- 
lem is to replace the interaction potential by a delta-like 



pscudo potential that supports only a single bound state 
and can be adjusted to have the same s-wave scattering 
length as the full interaction potential [l||. In this case 
the problem can be tackled, e.g., by choosing a basis of 
multi-band Wannier functions of the optical lattice 
However, effects of higher partial waves or of an energy 
dependence of the scattering length, that can easily be- 
come important if multiple Bloch bands are occupied [13] , 
are ignored by this approach. Here, the problem is ap- 
proached by using the stationary solutions of the lattice 
potential obtained by a procedure presented in (l8| . For 
this, the Hamiltonian is first separated into relative (rel.) 
motion and center-of-mass (cm.) motion. The different 
length scales are covered by expanding the rel. and cm. 
wave functions in spherical harmonics and a flexible basis 
of B splines for the radial part [l8|. In a configuration- 
interaction procedure the eigenfunctions of the rel. and 
cm. part of the full lattice Hamiltonian are used to de- 
termine the full eigenfunctions. These eigenfunctions are 
subsequently used as a basis for the propagation of the 
time-dependent wavefunction. 

The paper is organized as follows. In Sec |H] the sta- 
tionary Hamiltonian of the system is presented. In order 
to understand the numerical approa ch, the basis func- 
tions obtained by the procedure in [18[ are shortly in- 
troduced while the interested reader should consult [lH 
for a more detailed description. In Sec. Mil the time- 
propagation method is described. Afterwards in Sec IIVI 
the results of the time propagation are validated by a 
comparison to problems that possess an analytical solu- 
tion. Finally, in Sec. [V] the numerical method is used 
to analyse a system of 6 Li- 7 Li in a three-well OL. The 
experimentally relevant cases of an acceleration of the 
lattice, i.e., a linear perturbation, and of an additional 
harmonic confinement are considered. 



II. STATIONARY HAMILTONIAN AND ITS 
EIGENSOLUTIONS 



The full Hamiltonian 



H(t) =H + W{t) 



(1) 
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consists of a time-dependent part W(t) (specified below) 
and a stationary part 

(2) 

for two particles i = 1,2 with mass to^ interacting via 
the potential Vint- In the case of ultracold atoms the 
isotropic interaction potential Vint(^i — ? 2 ) = Vi n t(|fi — 
r*2 |) is described by an often only numerically given Born- 
Oppcnhcimer potential. The trapping potential 

V$= E V^ S m 2 (k uUi ) (3) 

u—x.y.z 

is that of an OL formed by three counter-propagating 
laser beams with wave vector k u in u direction (u = 

x, y, z). The lattice depth Vu' 1 '* is proportional to the laser 
intensity in direction u and the polarizability of particle 
i. 

The infinite lattice potential Vj a t is reduced to a poten- 
tial Viat with finite extension by expanding Vq a t to some 
arbitrary order into a Taylor series in all three directions 
[see Fig. Q] for the example of a 22nd order expansion of 
V x sin 2 (k x x)]. For practical purposes only expansions of 
order 2(2n + 1) are important, which lead to lattice po- 
tentials Viat with Viati?) — > oo for \r\ — > oo. For other 
expansions the potential is unbound from below leading 
to the appearance of non-physical states. 



The symmetry operations of D 2 h arc 

S = {E, C 2 (x), C 2 (y), C 2 (z), a(xy), cr(xz), <r(yz), i} , 

(4) 

where E is the identity, C n {u) is the rotation about — 
around the u axis (u = x, y, z), a{u\u 2 ) the reflection on 
the (iti, u 2 ) plane and i the inversion (i. e. point reflection 
at the origin). 

Since the interaction potential VS nt is invariant under 
any operation in S also the full unperturbed Hamilto- 
nian Hq belongs to the D 2 h point group if the symmetry 
operations are performed on both coordinates f± and f 2 
simultaneously. 

The group D 2 h possesses eight irreducible representa- 
tions r CT with 

a e {A g , Big, B 2g , B 3g , A u , Bi u , B 2u , B 3u } . (5) 

The characters of these irreducible representations are 
listed in Table Q] 
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FIG. 1: (color online) The 22nd-order expansion V\ a t(x,y = 
0, z = 0) of the lattice potential Viat (x, y = 0, z = 0) in IE- 
direction (solid line). Lengths are given in units of the lattice 
spacing d — n/k x . A linear perturbation as it appears, e.g., 
for an acceleration of the lattice in x direction leads to an 
inclination of the lattice sketched by the red dashed line, while 
an additional harmonic confinement raises the left and right 
lattice site (blue dotted line). 

The trapping potential Viat of an OL (and also Mat) has 
orthorhombic symmetry, which is characterized by the 
point group D 2 h- By adapting the basis functions to this 
symmetry, the eigenfunctions and the time-dependent 
wave function can be determined more efficiently. The 
symmetry of the problem is discussed in depth in [l8| . 
Here, only the essential points are repeated. 



TABLE I: Character table of the D^h point group. 

In order to find the eigensolutions of Hq the system is 
split into rel. and cm. coordinates, 



p = n — r 2 



R = 



m,\r\ + m 2 ?"2 



(0) 



mi + m 2 

With this separation, the Hamiltonian is written as 
H (R, p) = H c . m . (R) + H ml (p) + H coupl (R, p) , (7) 

where H c .m.i H ie \_, and H conp \ (i?, p) still have D 2 h~ 
symmetry [18| . 

The eigenfunctions of rel. and cm. are described in 
spherical coordinates and expanded in a basis of B splines 
B a and spherical harmonics Y" 1 . Since the symmetry 
operations of D 2 h commute with the Hamiltonian, the 
eigenfunctions can be chosen such that their symmetry 
properties correspond to some irreducible representation 
T a of D 2 h- In the following the rel. (cm.) eigenfunctions 



are denoted as <^ ct) (p) [^f'(R)} with j = 1, 2, 3, . . . 

In a configuration-interaction procedure products of 
eigensolutions of H CA11 , and /f rc i., i.e. configurations, 
are used to diagonalize the full Hamiltonian Hq. Be- 
cause all irreducible representations of D 2 h are one di- 
mensional, the direct product of two irreducible repre- 
sentations r K ® T\ is again an irreducible representation 



0), 



3 



T a that can be determined from the product table [TT1 
Hence, each configuration ^ (i?)<^- (/?) has the sym- 
metry properties of the related irreducible representation 
T a = T K <X> T\. The full solutions of a given symmetry a 
has the form of a superposition 



E 

{k,A}Gct ij 



E C &' A N K) 



(8) 



where {k, A} G a should indicate that the summation 
is performed over irreducible representations that fulfill 

r K ®r A = iv 

When considering identical bosonic (fcrmionic) par- 
ticles the rel. wavefunction has to be symmetric (an- 
tisymmetric) under inversion, i. e. only basis functions 
of rel. motion with A G {A g , B lg , B 2g , B 3g } (A G 
{A u , Bi u , B 2u , Bz u }) are used to form configurations. 



The wavefunctions, i. e. the coefficients C^'^ in 
Eq. ((SJ, arc finally determined by solving the eigenvalue 
problem 



= E {a) <J> (<t) 



(9) 



of Hq in the configuration basis. 
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TABLE II: Product table of irreducible representations of the 
D 2 h point group. 



for the evolution of the time-dependent coefficients 
B K j(t), which is governed by the matrix elements 



,(re,o-) _ 



<I> 



(«) 



W(t) 



) of the perturbation. Con- 



sidering two eigenstates 

*L T) )= E E c 

{k,A}£t ij 

*P)= E E c 

{^.,^}Gcr kl 



(«,A) 
ij 



kl 



A") 



(13) 



as specified in Eq. (|8]h the matrix elements of a pertur- 
bation are 



={*ti ) \w{t)\^) 



E E E E(c 

{k,A}Et ij {fi,u}ery kl 



(k,A) 



C 



kl 



(14) 



In general, the perturbation W'(t) can be expanded in 
a time-dependent Taylor series of its spacial coordinates 

w(t) = J2 E E /tf''(«ff. 

nm u—x,y,z u'—x,y,z 

where p u (i? u ) is the component of the rel. (cm.) motion 
in u direction (u = x, y, z). 

At the present stage perturbations in x direction of the 
general form 



W(t) =/oi(t) p x + /io(i) Rx + fu(t) R X R X 
+ f02(t) pI + f 20 (t) R 2 X 



(15) 



III. TIME-DEPENDENT EVOLUTION 

The Schrodinger equation of the time-dependent evo- 
lution 

&o+ft®)mt))=i*§ i m)) (10) 

with |* (t = 0)) = |* ) 

is solved in the basis {$i Cr ' > } of eigenfunctions of Hq of 
Eq. ©, 



|*(t)> = 5>o<(*) 



(<?) 



(11) 



Plugging Eq. (fTTj) into Eq. (|10p and multiplying from 



the left by ( $ 



leads to the equation 



(a) 



are implemented. In principle, the method can be easily 
extended to allow for perturbations in other directions 
and of higher orders. 

In order to illustrate how the perturbation matrix 
is computed, the case of a linear perturbation W = 
fio(t)R x is discussed in more detail. This perturba- 
tion does not couple the orthonormal rel. basis functions 



Thus, the summations in Eq. (|14[) reduce to 



vt^=ho(t) E E E E( c 



{k,A}€t ij {n,\}ea k 



M K) \RMk ] 



,(k,A) 

'ij 



C 



(M,A) 
kj 



,(«) 



R x 



(m) 



A- 



(16) 



is consid- 



In the following the term {^>\ 
ered for the exemplary case of n = A g . In this case the 
wave function ^ (i?) is totally symmetric (see Table|T]). 
Hence, ^/\^\r) needs to be anti-symmetric in x direc- 
(12) tion and symmetric otherwise, which is fulfilled solely for 
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jj, = i?3 U . In all other cases the integral vanishes. The 
according cm. basis functions are represented as 



^ 9) (i?,e,$) 



N a N t I 

E E E 

a=l /=0,{2} m=0 5 {2} 



(A.) i Ug) 



E E E 

a=l J=l,{2} m=l,{2} 



(B 3 „ 



1 J? 



Im ' 



(17) 



where B a are B splines, 



r z m (e,$) ± yr m (e,$) 



are a sum of spherical harmonics for m 7^ 0, and = 
Y" ; (O,$) (see [HI for details). The numbers in curly 
brackets below the sums indicate the summation step. 
N a and Ni are variable values of the number of B splines 
and the maximal angular momentum, respectively. With 
R T = R sin 6 cos $ one finds 



M k) \r x \^ 



(/<> 



1 v 

E E 

Z=0,{2} m=0,{2} Z' = l,{2} m'=l,{2} aa' 



E E E 



f c (As) V 
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'0 Jo 



Using the identities {Y[ 



) sin 9 cos ■ 
(18) 

(-l) m Yp", sin 8 cos $ = 



[yf 1 (e,$)-r 1 1 (e,$)],and 



p7T p2t7T 

/ sinede / d$r™ i (e,$)F™ 2 (e,$)r™ 3 (e,$) 

Jo Jo 



(2fi + l)(2Z 3 + l)(2/ 3 + l) /Zi Z 2 k 
Am lo 



Zi Z2 Z3 

TTil ?Ti2 7T13 

(19) 

the integral over the angles in Eq. (jT5)l can be efficiently 

^1 ^2 ^3 



computed in terms of Wigner 3i-symbols 

\mi rai m 3 

The other types of perturbations in Eq. (fT5|) are treated 
in an analogous way. 

Since the system is six dimensional the analy- 
sis in terms of the full time-dependent wavefunc- 
tion is nontrivial. However, equipped with the ma- 



trix elements of all perturbations, ($ 



RxPx 



(/<> 



(«) 



Pa 



(«) 



Px 



R, 



00 



fe 

, and 



fc / one can easily determine the expec- 
tation values of some of the most important observables. 
For example, the squared mean particle distance in x di- 
rection is given as 

(pi) =<*(*) I /£|*(f)) 



(<t) 



~2 



4' 



Likewise, one can determine the mean particle position 
or the uncertainty of the position in x direction. 



IV. COMPARISON WITH ANALYTICAL 
RESULTS 



In order to validate the numerical procedure a compar- 
ison with analytical results is necessary, which are avail- 
able for the harmonic approximation of the OL poten- 
tial. In the case of two identical particles of mass m in 
a harmonic trap the system decouples into rel. and cm. 
motion with Hamiltonian 

^° = Wi + \ Muj2r2 + Y + \ muj2p2 + v ^ ■ ( 2 °) 



Here, M = 2m, /1 = m/2, P is the momentum of cm. 
and p the momentum of rel. motion. In the following, 
a linear perturbation W(t) = f{i)R x and a quadratic 
perturbation W(t) = f(t)R%., i.e. a time-dependent ac- 
celeration and a variation of the trapping frequency, are 
considered. Since the cm. part of Hq decouples into x, y, 
and z direction, only the cm. harmonic oscillator in x 
direction with Hamiltonian 



P 2 



1 



P 2 



1 Rl 



is affected by the perturbations, where Ah = y/h/ (McJ) 
is the harmonic oscillator length. The comparisons are 
performed for expectation values of the position 



X(t) = (R x 



*(t) R x *(t) 



and the mean deviation from X 



v(t) = \ (R 2 X ) - (R x 



A. Periodic driving 



(22) 



(23) 



For the case of a periodically driven harmonic oscillator 
with driving strength C s hake and frequency wo, 

R 

W 1 (t) = Cshake cos (w t) -p- , (24) 
there exists an analytic solution fl9j . 

i> n {X,t) = e i ^<j )n {X-m), (25) 

where ip{x, t) is a phase, which vanishes for t = 0, <p n is 
the nth harmonic oscillator eigenstate of -ffho; and 



/- / i\ ^ho Cshake / . \ 

= 1 2T^cos(w t). 

1 — 0Jq/0J z 



(26) 
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In order to conform with the initial condition 

ip n (x,o) = <p n (x - m) 



(27) 



the trap is shifted at t = to £(0) by instantly adding a 
constant linear perturbation 



W 9 = -hu>a 



1 



hake 



l-^ /^A ho - (28) 



From the analytic solution one obtains straightfor- 
wardly 

X(t) = -A ho C BhakD [1 - cos(wot)] , a(i) = %. (29) 

v2 

In Fig. [5] a comparison of a numerical calculation of 
X(t) to the result in Eq. ([29|) shows very good agree- 
ment with deviations on the order of 10~ 10 . A similar 
accuracy is obtained for the value of <j(t). The devia- 
tions are due to the finiteness of the basis which, in the 
shown calculation, only includes basis functions with an 
cigenencrgy below the chosen cutoff of 20hui. The en- 
ergy cutoff can be adapted to reach higher accuracies, if 
needed. 
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Thus, assuming perfect adiabaticity, the width of the 
wave function behaves like 



<y{t) = A ho / V2(l + 2C harm wt) 



(31) 



In Fig. [3] a comparison to the numerical calculations 
shows good agreement to this result with an error of 
about 5 x 10 -5 for Charm = 0.002, which is due to nona- 
diabatic effects. For example, reducing the speed of the 
perturbation by setting Charm = 0.001 reduced the error 
to about 2 x 10~ 5 . 
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FIG. 3: (color online) Comparison of analytical (blue solid) 
and numerical (black dashed) results for cr(t) [see Eqs. (|23|) 
and ([31) ] for Charm = 0.002. The error A error = \a — <r num | is 
shown in the inset. The relatively large error in comparison 
to the results shown in Fig. [2] is due to nonadiabatic effects. 
These effects get smaller for larger t since the change of A^ a (t) 
is reduced [see Eq. (|30[) ]. For cut > 5000, however, the incom- 
pleteness of the basis used for the numerical calculations (only 
states with energies below E = 20 hu) are included) leads fi- 
nally to an increase of A crror . 



FIG. 2: (color online) Comparison of analytical (blue solid) 
and numerical (black dashed) results for X(t) [see Eqs. (|22p 
and (|29|) ] for C S h a ke = 0.5. The difference of the results is 
below 10 -10 and therefore invisible. The width of the wave 
function a(t) = Aho/v2 is numerically reproduced with the 
same level of accuracy. 



B. Adiabatic deepening 

The mean width of the wavefunction a for an har- 
monic oscillator with oscillator length Ah Q is given as 
Aho/V% [see Eq. (|2"9"|) ]. Considering a time dependent 
perturbation W(t) = C\ laim HujR^,/A^ 10 ujt, the full poten- 

1 Ft 2 

tial is given as ^hjjj-^-{l+2Ch arm ujt). If the perturbation 
happens sufficiently slowly, the wave function will always 
remains in an cigenstatc of a harmonic oscillator with a 
trap length 



A ho (t) = A ho (t = 0)/y/l + 2C haxm ut. 



(30) 



V. EXAMPLE CALCULATIONS FOR 6 LI- 7 LI 

In the following a system of two distinguishable atoms, 
6 Li and Li is considered. The interaction potential Vint 
is given by the Born-Oppcnhcimcr potential for scat- 
tering of spin-polarized lithium. The atoms are con- 
fined in a three-site lattice potential T4at> which is re- 
alized by a 22nd order expansion of Hat m Eq. © 
in x direction (see Fig. [1]) and a harmonic approxima- 
tion in y and z direction. The chosen wave vectors 
k x = k y = k z — 27r/(1000nm) lead to a lattice spacing of 
d = 500 nm = 9450 a.u. A lattice depth in x direction of 
V x = 1.36huji, where uq is the frequency of the harmonic 
approximation of the lattice for atom 1 ( 6 Li), results in 
the relatively small hopping energies J\ = 2.1 x \Q~ A ftwi 
of atom 1 and Ji = 1.5 X \Q~ A ftwi of atom 2 in the corre- 
sponding Hubbard model for the infinite lattice. Hence, 
even for the relatively small s-wave scattering length of 
41 a.u. of 6 Li- 7 Li a correlated Mott-like state is formed, 
i. e., the atoms do not occupy the same lattice site in the 
ground state [2l|. Since no unit filling of the lattice is 



considered, the atoms are nevertheless mobile in x direc- 
tion. This enables the observation of a correlated motion 
of the distinguishable atoms. The lattice depths in y and 
z direction are given as V y = V z w 8V X such that for 
low-lying states motion in these directions is frozen out. 

Despite the reduction to only three lattice sites, the 
considered system exhibits the basic mechanisms of hop- 
ping and onsite interaction of atoms in an OL. Similar 
systems of only a few lattice sites appear also experimen- 
tally in superlattices fic| . 



A. Linear perturbation 

First, the system is adiabatically inclined by a pertur- 
bation of the type W(t) = AtR x [see Fig. [4] (a)]. Experi- 
mentally this could, e.g., be realized by slowly increasing 
the acceleration of the lattice in x direction. The system 
starts in the ground state where the atoms spread sym- 
metrically over the lattice. As a consequence, the mean 
atom position is exactly in the middle of the triple- well 
potential, i.e. at x/d = 0. Due to their repulsion the 
atoms never occupy the same lattice site. In this case 
their mean distance \J (p x ) is approximately d. The cor- 
responding probability density along the x axis is shown 
in the left graph of Fig. @] (b). 

Upon inclining, the system stays in the state of mini- 
mal energy, i. e. the heavier 7 Li atom slowly moves into 
the lower left lattice site (i.e. x 2 = (x 2 ) approaches —d) 
while the lighter 6 Li atom moves to the central site (i. e. 
x\ = (.ii) approaches zero), where it avoids an energy 
gain due to the interatomic repulsion. With much smaller 
probability the same process with exchanged 6 Li and 7 Li 
appears [see right graph of Fig.|4](b)]. During the process 
the mean distance is unchanged while the uncertainty of 
the position \J (Jxi — Xi) 2 ) of atom i (i = 1,2) decreases 
[see Fig. |4] (a) and (b)]. Stopping at a final inclination 
that results in an energy difference of 0.09fiw between 
neighboring wells, the atoms are well separated. For a 
further inclination both 6 Li and 7 Li would move to the 
left well. 

Starting from a system of separated atoms, one can 
induce a collision process. To this end the linear pertur- 
bation, i. e. the acceleration, is suddenly switched off. As 
shown Fig. 2] (c) in this case the heavier atom tunnels 
back and forth between the left and the right well. Due 
to the small initial population of the state where 6 Li is 
in the left well and 'Li in the central well, also 6 Li tun- 
nels back and forth and x\ oscillates slightly around zero. 
Owing to the mass difference both tunneling processes of 
6 Li and 7 Li happen with different frequencies. Due to the 
repulsion during the tunneling process the atoms do still 
not occupy the same lattice site which is obvious from 
the unchanged particle distance. 

While a weak adiabatic inclination can be easily de- 
scribed also within the standard Hubbard model, a fast 
inclination couples states of different Bloch bands 0] . In 
Fig. [3] the behavior for a stronger and faster inclination 
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FIG. 4: (color online) Mean particle position Xi = (xi) of 6 Li 
(thick lighter red line) and of 7 Li (thick darker blue line) and 
mean distance *J (p|) (grey dashed line). The corresponding 
lighter red, and darker blue shading illustrates the uncertainty 
of the position Xi ± yj (Jxt — Xi) 2 ) of 6 Li and 7 Li, respectively. 
Time is given in units of the hopping time H/Ji of 6 Li. (a) 
Time dependent behavior for a linear inclination with a final 
perturbation W = 66,hR x /d = 0.09 hoiR x /d. (b) Probabil- 
ity density |^(xi, X2)\ for yx = y 2 = Zi = 22 = of the 
initial state (left) and the final state (right). Initially there is 
an almost equal probability of finding 6 Li in the central well 
and 7 Li in the outer wells and vice versa. After the linear in- 
clination 7 Li is predominantly situated in the left well and 6 Li 
in the central well. The situation with exchanged 6 Li and 7 Li 
the wavefunction has a small but non-vanishing probability, 
(c) Free evolution of the system with the initial state being 
the final state of the process in (a). 
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than the one in Fig. [4] (a) is presented. In this case the 
behavior is harder to predict. For example, it is unclear 
whether cither first the heavier atom or the lighter atom 
moves to the left lattice site. Although one could ex- 
pect that the lighter atom with its larger tunneling rate 
is more mobile and will move first, indeed the heavier 
atom tunnels first to the left well. During the fast in- 
clination also states with two atoms at the same lattice 
site are occupied, w hich is reflected by a reduction of 
the mean distance y/ (p%.). The occupation probability 
of states above the first Bloch band is high (see top of 
Fig- 13: an d thus the behavior cannot be described within 
a single-band approximation of the Hubbard model. 
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FIG. 5: (color online) Time dependent behavior for a linear 
inclination with a final perturbation W = 4400Ji_R :E /d = 
0.6 hw\R x jd. Top: Total occupation probability of states 
above the first Bloch band. Bottom: Legend as in Fig. [4] 



B. Harmonic perturbation 

In experiments optical lattices are not infinite but 
the atoms are normally confined by an additional weak 
harmonic potential. In the following the effect of the 
sudden activation of such a harmonic potential W = 
A(x\ + x^/d 2 is studied. This perturbation does not 
break the symmetry of the potential and the mean po- 
sition of the atoms remains at x/d = 0. However, as 
one can see in Fig. [5] (a) for a certain strength of the 
harmonic perturbation the system oscillates between un- 
bound states (V\pI) ~ d) and repulsively bound states 
(V (Px) ~ 0.5c£) [13] that are in resonance. These oscil- 
lations are also visible in the uncertainty of the atoms' 
positions. For an increased harmonic perturbation no re- 
pulsively bound state is in resonance with the unbound 
state. Hence, as shown in Fig. HJb) the atoms oscillate 
predominantly between delocalized states and states lo- 



calized at the central lattice site. Since the atoms re- 
pel each other, the oscillations are exactly opposing each 
other. The off-resonant coupling to the bound state l eads 
to small and fast oscillations of the mean distance y/ (p 2 ,) 
between 0.8d and l.Od. 
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FIG. 6: (color online) Time-dependent behavior for the sud- 
den turn-on of an additional harmonic confinement, (a) For 
W = 82 Ji(xl + xf)/d 2 = 0.01&Ji(f? + il)/d 2 oscillations 
between bound and unbound states appear, (b) For stronger 
confinement W = 164 J L (z? +x 2 2 )/d 2 = 0.02 fiux{xl + xl)/d 2 
the bound-state occupation is much weaker, however the par- 
ticles tunnel alternating between the central and outer wells. 
Legend as in Fig. [4] 



VI. CONCLUSION AND OUTLOOK 

A theoretical approach for the full non-pcrturbative 
time-dependent description of two interacting particles 
in an optical lattice was introduced. A comparison 
with analytical results shows the possibility to perform 
high-precision analyses. Example calculations for 6 Li- 7 Li 
in a three-well optical lattice where performed, demon- 
strating the possibility to analyze this complex six- 
dimensional system in terms of several expectation val- 
ues. It was shown how the atoms are separated by a 
slowly increasing acceleration of the system and how the 
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system reacts upon suddenly stopping the acceleration. 
It was also demonstrated that a fast acceleration of the 
lattice leads to a strong occupation of states above the 
first Bloch band, which marks the break down of the 
usually adopted single-band Hubbard models. As finally 
shown, a weak harmonic perturbation can have an impor- 
tant impact if the system encounters a resonance between 
bound and unbound states. 

The use of a spectral method, i. e. expanding the time- 
dependent wavefunction in a basis of eigenfunctions of 



some underlying Hamiltonian, offers a large degree of 
flexibility. For example, by modifying the underlying 
Hamiltonian the here-presented system of two neutral 
atoms can be easily generalized to other particles, such 
as ions, or dipoles. Also the external potential is flexible 
enough to describe a large class of systems like quantum 
dots or one- and two-dimensional optical traps. In the 
future we intend to analyse and develop with the pre- 
sented procedure schemes for the fast and high-fidelity 
manipulation of small quantum systems. 
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